#ifndef EQUATIONS_H
#define EQUATIONS_H

#include "signal.h"

// -------------------------------------------
// C2 = C2(T,W)
// INPUTS:
// W: nm/s
// T: K
// OUTPUTS:
// C2: nm^2/s

class eq_C2:public MathFunction{
 private:
  // parameters
  int LatticeSize;
  double XYmax;
 public:
  eq_C2(int i_LatticeSize = 100,double i_XYmax=8000);
  virtual ~eq_C2();
  virtual double value(const double *x);
  virtual void gradient(const double *x,double *grad);
  virtual bool parameterize(const double *para);
  virtual bool parameterize(int para_ID,double para_val);
};
// ------------------------------------------
// sigma^2 = sigma^2(W)
// INPUTS:
// W: nm/s
// OUTPUTS:
// sigma^2: nm^2/s^2

class eq_sigma2:public MathFunction{
 private:
  int LatticeSize;
  double XYmax;
 public:
  //  eq_sigma2(int i_LatticeSize = 100);
  eq_sigma2(int i_LatticeSize=100, double i_XYmax=8000);
  virtual ~eq_sigma2();
  virtual double value(const double *x);
  virtual void gradient(const double *x,double *grad);
  virtual bool parameterize(const double *para);
  virtual bool parameterize(int para_ID, double para_val);
};
// ------------------------------------------
// cov_z_p = cov_z_p(W) = cov_z_p(C2(W),sigma2(W))
//class eq_cov_z_p:public MathFunction{
// private:
//  // parameters
//  int p;
//  int m;
//  int n;
//  double dt;
//  double cov_z_m;
//  double XYmax;
//  int    LatticeSize;
//  // Auxillary variables
//  int index;
//  eq_C2 f_c2;
//  eq_sigma2 f_sigma2;
// public:
//  eq_cov_z_p(int i_p=1,int i_m=5,int i_n=0,
//	     double i_XYmax = 8000,int i_LatticeSize = 100,
//	     double i_dt=0,double i_cov_z_m=0);
//  virtual ~eq_cov_z_p();
//  virtual double value(const double *x);
//  virtual void gradient(const double *x,double *grad);
//  virtual bool parameterize(int para_ID,double para_val);
//};
// ----------------------------------------
// z2_p = z2_p(W,A) = z2_p(C2(W),sigma2(W),W_pmn(A))
class eq_z2_p:public MathFunction{
 private:
  // parameters
  int p;
  int m;
  int n;
  double dt;
  double cov_z_m;
  double mean_z_m;
  double Freq;
  double XYmax;
  int    LatticeSize;
  // auxillary variables
  int index;
  eq_C2 f_c2;
  eq_sigma2 f_sigma2;
 public:
  eq_z2_p(int i_p=1,int i_m=5,int i_n=0,
	  double i_dt=100,double i_cov_z_m =0,double i_mean_z_m=0,
	  double i_Freq=5, double i_XYmax=8000,int i_LatticeSize=100);
  virtual ~eq_z2_p();
  virtual double value(const double *x);
  virtual void gradient(const double *x,double *grad);
  virtual bool parameterize(int para_ID,double para_val);
  virtual bool parameterize(const double* para);

  // Additional functions
  double eq_W_pmn(double W,double A);
  double eq_dW_pmn_dA();
};
// --------------------------------------------
// cov_z_p = cov_z_p2(C2(W),sigma2(W))
//class eq_cov_z_p2:public MathFunction{
// private:
//  // parameters
//  int p;
//  int m;
//  int n;
//  double dt;
//  double cov_z_m;
//  // Constant
//  static const double XYmax;
// public:
//  eq_cov_z_p2();
//  virtual ~eq_cov_z_p2();
//  virtual double value(const double *x);
//  virtual void gradient(const double *x,double *grad);
//  virtual bool parameterize(int para_ID,double para_val);
//};
// --------------------------------------------
// r2_f = r2_f(W,A)
class eq_r2_f:public MathFunction{
 private:
  // Parameters
  int mode;
  double XYmax;
  double Freq;
  int LatticeSize;
  double dt;
  double *cov_z_m[4];
  double *mean_z_m[4];
  
  // Auxillary variables
  eq_z2_p z2_p;
 public:
  eq_r2_f(double i_dt = 100.0,int i_mode = 50,double i_XYmax = 8000.0,
	  int i_LatticeSize = 100,double i_Freq=5.0);
  virtual ~eq_r2_f();
  virtual double value(const double *x);
  virtual void gradient(const double *x, double *grad);
  virtual bool parameterize(int para_ID, double para_val);
  virtual bool parameterize(int para_ID,double *para_val);
};
// ------------------------------------------
// m2_f = m2_f(W,A)
class eq_m2_f:public MathFunction{
 private:
  // Parameters
  int mode;
  double XYmax;
  double Freq;
  int LatticeSize;
  double dt;
  double *cov_z_m[4];
  double *mean_z_m[4];
  double *K[4];
  
  // Auxillary variables
  eq_z2_p z2_p;
  void CalculateK();
 public:
  eq_m2_f(double i_dt = 100,int i_mode = 50,double i_XYmax = 8000,
	  int i_LatticeSize = 100, double i_Freq =5);
  virtual ~eq_m2_f();
  virtual double value(const double *x);
  virtual void gradient(const double *x,double *grad);
  //  virtual bool parameterize(const double *para);
  virtual bool parameterize(int para_ID,double para_val);
  virtual bool parameterize(int para_ID,double* para_val);
};
// ------------------------------------------
// f = f(W,A) = q_r*(r2_set-r2_f)^2+q_m*(m2_set-m2_f)^2
class eq_f:public MathFunction{
 private:
  double r2_set;
  double q_r;
  double m2_set;
  double q_m;
  int    mode;
  int    LatticeSize;
  double XYmax;
  double dt;
  double Freq;

  double *cov_z_m[4];
  double *mean_z_m[4];

  eq_r2_f f_r2_f;
  eq_m2_f f_m2_f;
 public:
  eq_f(double i_r2_set = 900,double i_q_r=1.0,
       double i_m2_set = 36,double i_q_m=1.0,
       double i_dt = 100, int i_mode=50,int i_LatticeSize = 100,double i_XYmax = 8000,double i_Freq = 5);
  virtual ~eq_f();
  virtual double value(const double *x);
  virtual void gradient(const double *x,double *grad);
  virtual bool parameterize(int para_ID,double para_val);
  virtual bool parameterize(int para_ID,double *para_val);
};

#endif
